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The phaseless auxiliary-field quantum Monte Carlo (AF QMC) method [S. Zhang and H. Krakauer, Phys. 
Rev. Lett. 90, 136401 (2003)] is used to carry out a systematic study of the dissociation and ionization energies 
of second-row group 3A-7A atoms and dimers, Al, Si, P, S, CI. In addition, the P 2 dimer is compared to the 
third-row As 2 dimer, which is also triply-bonded. This method projects the many-body ground state by means 
of importance-sampled random walks in the space of Slater determinants. The Monte Carlo phase problem, 
due to the electron-electron Coulomb interaction, is controlled via the phaseless approximation, with a trial 
wave function Y&t). As in previous calculations, a mean-field single Slater determinant is used as |^t). The 
method is formulated in the Hilbert space defined by any chosen one-particle basis. The present calculations use 
a planewave basis under periodic boundary conditions with norm-conserving pseudopotentials. Computational 
details of the planewave AF QMC method are presented. The isolated systems chosen here allow a system- 
atic study of the various algorithmic issues. We show the accuracy of the planewave method and discuss its 
convergence with respect to parameters such as the supercell size and planewave cutoff. The use of standard 
norm-conserving pseudopotentials in the many-body AF QMC framework is examined. 

PACS numbers: 71.15.-m,02.70.Ss,31.25.-v,31.15.Ar 



I. INTRODUCTION 

Achieving accurate solutions of the electronic many-body 
Schrodinger equation is a challenging problem for calcula- 
tions of the properties of real materials. For many systems, 
density functional theory (DFT), in a variety of approxima- 
tions, has been applied with great success. In DFT, the many- 
body interactions are replaced by a single particle interacting 
with the mean-field generated by the other particles, similar in 
spirit to the Hartree-Fock (HF) method. Unfortunately, these 
methods have well-known limitations and often fail at describ- 
ing the properties of materials with large electron-electron 
correlation. 

A more accurate approach is the quantum Monte Carlo 
(QMC) metho d 1,2 ' 3 ' 4 , which has been shown to be among 
the most effective methods for many-electron problems. Un- 
like other correlated methods, QMC calculation times scale 
as a low power of the system size^. The fixed-node diffusion 
Monte Carlo (DMC) approach, which samples the many-body 
wave function in real space, has been the most widely used 
QMC method in electronic structure calculations^ 3 -. 

The recently developed phaseless auxiliary-field quantum 
Monte Carlo (AF QMC) method 4 - is an alternative and com- 
plementary QMC approach, which samples the many-body 
wave function in the space of Slater determinants. This 
method has several attractive features. The fermionic anti- 
symmetry of the wave function is automatically accounted for, 
since it is sampled by Slater determinants. This provides a 
different route to controlling the sign problem^^ from fixed- 
node DMC, which has shown promise in reducing the depen- 
dence of the systematic errors on the trial wave functio n 9,10 ' 11 . 
The orbitals in the Slater determinants are expressed in terms 
of a chosen single-particle basis {e.g., planewaves, Gaussians, 
etc.), so AF QMC shares much of the same computational ma- 
chinery with DFT and other independent-particle type meth- 
ods. AF QMC can thus straightforwardly incorporate many of 



the methodological advances from mean-field methods (such 
as pseudopotential and fast Fourier transforms) while system- 
atically improving on mean-field accuracy. 

Using a planewave basis, tests of the phaseless AF QMC 
method for a few simple atoms and molecules 4 ^ 2 - as well as for 
the more correlated TiO and MnO molecules^ yielded excel- 
lent results. More systematic applications of the phaseless AF 
QMC method to atoms and molecules have been carried out 
using Gaussian basis sets. These include all-electron calcula- 
tions for first-row systems'^ and effective-core potential cal- 
culations in post-<i group elements^. The results also showed 
excellent agreement with near-exact quantum chemistry re- 
sults and/or experiment. 

The planewave AF QMC method is well-adaped for corre- 
lated calculations of extended bulk systems, where planewave 
based methods have been the standard choice in traditional 
electronic structure calculations. It is therefore important to 
systematically study its algorithmic issues and to characterize 
its performance. In this paper, we use the planewave phase- 
less AF QMC method to carry out a systematic study of the 
dissociation and ionization energies of second-row atoms and 
dimers in Group 3A-7A, namely Al, Si, P, S, CI. The inter- 
esting case of the triply-bonded P2 dimer is also compared to 
the third-row AS2 dimer. The principal goal of this study is 
to further benchmark the AF QMC method across more sys- 
tems and across different basis sets and to compare the re- 
sults with those from other methods and experiment. While 
the use of localized basis sets, such as Gaussians, is generally 
more efficient for isolated atoms and molecules, it is straight- 
forward to apply planewave methods using periodic boundary 
conditions and large supercells. Planewave methods have sev- 
eral desirable features. A planewave basis provides an unbi- 
ased representation of the wave functions, since convergence 
to the infinite basis limit is controlled by a single parameter, 
the kinetic-energy cutoff E cu[ . Planewaves are algorithmically 
simple to implement, and operations with planewaves can be 
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made very efficient, using fast Fourier techniques as in DFT 
methods. To keep the planewave basis size tractable, pseu- 
dopotentials must be used to eliminate the highly localized 
core electron states and to produce relatively smooth valence 
wave functions. The present choice of isolated atomic and 
molecular systems permits direct comparisons with Gaussian- 
based AF QMC and with well-establised quantum chemistry 
all-electron and pseudopotential calculations. 

The remainder of this paper is organized as follows. In Sec- 
tion!]]] we describe computational details of the phaseless AF 
QMC method in the planewave-pseudopotential framework. 
In Section [Hi] we discuss calculation parameters such as su- 
percell (simulation cell) size, cutoff energy, and AF QMC 
time step size, and the algorithm's convergence behavior with 
respect to these parameters. Section [IV] presents the calcu- 
lated dissociation and ionization energies and comparisons 
with other theoretical results and with experiment. In Section 
[VJ we discuss systematic errors due to the use of the phaseless 
approximation and use of norm-conserving pseudopotentials. 
We then conclude with some general remarks in SectionlvTl 



II. PLANEWAVE AF QMC METHOD: COMPUTATIONAL 
DETAILS 

A. Hamiltonian 

It is convenient to express, within the Born Oppenheimer 
approximation, the electronic Hamiltonian in second quan- 
tized form in terms of a chosen orthonormal one-particle basis 
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where M is the number of basis functions, cj and a are 
the corresponding creation and annihilation operators, and the 
electron spins have been subsumed in the summations. H^' 
and Hyjy are the one- and two-body matrix elements, and 
V\\ is the classical Coulomb interaction of the point ions 13 . 
Atomic units are used throughout this paper. We use periodic 
boundary conditions and a planewave basis 
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where fl is the volume of the simulation cell and G is a recip- 
rocal lattice vector. As in planewave-based density functional 
calculations, the number of planewaves M in the basis is de- 
termined G 2 /2 < E cut , where E'cut is the cutoff kinetic energy. 

The one-body operators in the Hamiltonian include the ki- 
netic energy, 
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and nonlocal pseudopotential, which describes the electron- 



ion interaction 
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where VL(G — G') and Vnl(G, G') are the matrix elements of 
local and nonlocal parts of the pseudopotential, respectively. 
It is convenient to rewrite the local part of the pseudopotential 
and to define the the following quantities 
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t4i,L = \ Y MQ) [p(Q) + PHQ)] , (5b) 
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where N is the number of electrons, and the one-body density 
operator /5(Q) in this equation is given by 

p(Q) = E c g+q,a c g,a (Ecut - |G + Q| 2 /2) , (6) 

G,A 

where the step function ensures that (G + Q) lies within the 
planewave basis, and the summation over electron spins (A = 
1, 2) has been made explicit. 

The electron-electron interaction is given by 
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The primed summation indicates that the G; = G/. singu- 
lar term is excluded, due to charge neutrality. The first term 
in this equation is a constant due to the self-interaction of an 
electron with its periodic images. It depends only on the num- 
ber of electrons in the simulation cell and the Bravais lattice 
associated with the periodic boundary conditions. The stan- 
dard Ewald expression for £ is given by 14 
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where R is a direct lattice vector, and £ is independent of 
the Ewald constant k, which only controls the relative con- 
vergence rates of the direct and reciprocal space summations. 
For the discussion below, we rewrite the two-body contribu- 
tion in Eq. (0: 
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The third term in Eq. (O is a sum of diagonal one-body opera- 
tors arising from the anticommutation of the fermion creation 
and destruction operators. 

Finally, we can regroup the contributions to the Hamilto- 
nian into constant, one-body, and two-body parts, 
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It is convenient to express the two-body part as a quadratic 
form of one-body operators (this can always be done, and such 
forms are not unique). We use the identity p{— Q) = pl(Q) 
to write 



^ (2) = E 07T2- \m)p\Q) +/5 f (Q)p(Q)] 
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Defining Hermitian operators A(Q) and S(Q) as 
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the two-body contribution becomes a simple sum of quadratic 
operators, 
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where we have used the Q 
expression. 
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B. Ground state projection and the Hubbard-Stratonovich 
Transformation 



The ground state of -ffl^o) = -Ebl^o) is obtained by 
imaginary-time projection from a trial wave function | *S?t) 
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provided (^xl^o) 0- 1 R tne present calculations, |^x) is 
a single Slater determinant obtained from a mean-field calcu- 
lation. Expressing the imaginary-time projection in terms of 



the small discrete time step At facilitates the separation of the 
one- and two-body terms, using the short-time Trotter-Suzuki 
decompositio n 15 ' 16 
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The application of the one-body propagator e ( 1 / 2 ) Ari ^ < ' on 
a Slater determinant |0) simply yields another Slater determi- 
nant: |0') = e^( 1 / 2 ) A ' ri ^ <1) |0). The two-body propagator is 
expressed as an integral of one-body propagators, using the 
Hubbard-Stratonovich transformationiii^ 
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for any one-body operators { b\ }. Thus we have 
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where we introduce a vector er = {<7i}, whose dimensionality, 
dim(er), is the number of all possible Q-vectors satisfying 
Q = G G' for two arbitrary wave vectors G and G' in the 
planewave basis. The operators v = { >/— A7 bi} are given by 
the iA(Q) or i-B(Q) one-body operators, since all the Aj = 1. 

In the original formulations of the AF QMC metho d 19,20 , 
the many-dimensional integral over the auxiliary fields <x 
in Eq. ( fT8l is evaluated by standard Metropolis or heat- 
bath algorithms. We instead apply an importance-sampling 
transformatio n 4 ' 2 1 1 22 to turn the projection into a branching 
random walk in an overcomplete Slater determinant space. 
The importance sampling helps guide the random walks ac- 
cording to the projected overlap with the trial wave function. 
More importantly, it allows the imposition of a constraint to 
control the phase problem. 

A phase problem arises for a general repulsive two-body in- 
teraction, because the A^ cannot be made all negative. In other 
words, not all components of the operator v can be made real. 
(Although this is in principle possible by an overall shift to the 
potential 20 or by introducing many more auxiliary fields, they 
both cause large fluctuations 4 ^.) As the random proceeds, 
the projection in Eq. ( fT8l 
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by a complex v causes the orbitals in the Slater determinants 
|0) to become complex. For large imaginary projection times, 
the phase of each |0) becomes random, and the stochastic rep- 
resentation of the ground state |^o) becomes dominated by 
noise. This leads to the phase problem and the divergence 
of the fluctuations. The phase problem is of the same ori- 
gin as the sign problem that occurs when the one-body op- 
erators v are real, but is more severe because, instead of a 
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+ \<f>) and — |0) symmetry^ 2 !, there is now an infinite set 
{e l9 \4>) , 6 e [0, 2ir)}, among which the Monte Carlo sam- 
pling cannot distinguish. 

The phaseless AF QMC method 4 used in this paper con- 
trols the phase/sign problem in an approximate manner using 
a trial wave function. The method uses a complex importance 
function, the overlap (^>t\<I>), to construct phaseless random 
walkers, \4>) / (^t\4>)' which are invariant under a phase gauge 
transformation. The resulting two-dimensional diffusion pro- 
cess in the complex plane of the overlap (tyj<\<f)) is then ap- 
proximated as a diffusion process in one dimension. Addi- 
tional implementation details can be found in Refs. H101I241 
The phaseless constraint is different from the nodal condition 
imposed in fixed-node DMC, since the phaseless constraint 
confines the random walk in Slater determinant space accord- 
ing to its overlap with a trial wave function, which is a global 
property of \<j>). Thus, the phaseless approximation can be- 
have differently from the fixed node approximation in DMC. 

Finally, we describe the use of fast Fourier transform (FFT) 
with a planewave basis to efficiently implement the random 
walk projection given by Eqs. d 1 3a| >, ( fT4b . ( fT~8T >, and ( fT9l ). For 
example, 
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Terms in the series can be evaluated as an iterative FFT, since 
p(Q)|0) is just a convolution. For typical values of At, we 
find that n max ~ 4 accurately reproduces the propagator. 



Matrix elements of one-body terms in the local energy (and 
other similar estimators) can be expressed in terms of the one- 
body Green's function o 21 ' 22 
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The Green's function can be expressed in terms of the one- 
particle orbitals in the |\&t) an d 10) Slater determinants as 
follows. A general Slater determinant \<p) can be written as 
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where the 4 creates an electron in the orbital i 
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and j labels the one-particle orthogonal basis functions, which 
are planewaves in the present case. The are the elements 
of a M x N dimensional matrix <I>. Each column of the matrix 
4> represents a single-particle orbital expressed as a sum of 
planewaves. It is a well-known result that the overlap of two 
Slater determinants is given by the determinant of the overlap 
matrix of their one particle orbitals 
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Finally, it can be shown that the Green's function can be ex- 
pressed as 2 ^ 
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C. Ground-state mixed estimator 

The ground state energy Eq can then be obtained by the 
mixed estimator 

= (*Tia|*o) = Um mHe-ey T ) 

which is evaluated periodically from the ensemble of Slater 
determinants generated in the course of the random walks. 
In the phaseless AF QMC method, an importance sampling 
transformation 4 " leads to a stochastic representation of the 
ground-state wave function in the form of 
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Hamiltonian matrix elements of two-body terms in the 
mixed estimator are expressed in terms of the two-body 
Green's function, which can be written as products of one- 
body Green's functions using the Fermion anticommutation 
properties, 
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Rather than directly implementing Eq. ( f30b , it is more efficient 
to use fast Fourier transformations to take advantage of local- 
ity in real space. The computer time to calculate the mixed 
estimator then scales as N 2 M log(M), where N is the num- 
ber of electrons and M is the number of planewaves. 



This means the mixed estimate for the energy is given by 
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where the local energy is defined as 
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D. Trial wave function 

The trial wave function |*t) determines the systematic ac- 
curacy of our calculations, due to the use of the phaseless ap- 
proximation. Its quality also affects the statistical precision. 
We use a single Slater determinant as the trial wave func- 
tion, which is obtained either from a DFT calculation or HF 
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calculation. The DFT wave functions were generated self- 
consistently with ABINIT— , using a planewave basis and the 
local density approximation (LDA). The HF wave functions 
were obtained from an in-house planewave-based Hartree- 
Fock program. In both cases, identical setup is used in the 
independent-electron calculation as in the corresponding AF 
QMC calculations. 



E. Pseudopotentials 



TABLE I: Optimized _B cut and OPIUM pseudopotential parameters 
used in the calculations. Each angular component (I) of the pseu- 
dopotential has its own cutoff radius (r c ). 



Tc (units of ao) Reference 
Species E mt (Ha) 1 = 1 = 1 1 = 2 configuration 



Al 


7.50 










Si 


6.13 


2.20 


2.20 


2.50 


[Ne] 3s 2 3p 2 


P 


18.00 


1.75 


1.75 


2.50 


[Ne] 3s 2 3p 2 5 3d - 5 


s 


19.00 


1.75 


1.75 


1.75 


[Ne] 3s 2 3p 3 - 5 3d - 5 


CI 


18.00 


1.75 


1.75 


2.50 


[Ne] 3s 2 3p 4 5 3d - 5 


As 


18.00 


1.80 


1.80 


2.50 


[Ar] 4s 2 4p 2 5 4d° 



Norm-conserving pseudopotentials^. are used in the 
present calculations. Pseudopotentials are necessary to keep 
the basis size tractable by eliminating the highly-localized 
core states. Pseudopotential transferability is a source of po- 
tential errors, however, especially since the pesudopotentials 
used here are generated from independent-electron calcula- 
tions. Such pseudopotentials are quite routinely employed in 
QMC and other many -body calculations and have proved very 
useful. But their transferability is not nearly as extensively 
quantified and studied as in standard independent-electron 
calculations. Thus one of our goals here is to examine the use 
of such pseudopotentials in the many-body AF QMC frame- 
work. 

The pseudopotential has been adapted to take the 
Kleinman-Bylander (KB)2I form suitable for planewave cal- 
culations, 
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where tpi is the pseudoorbital for the Z-th angular momentum 
component. In our case, we use the neutral atomic reference 
state (with an LDA-type Hamiltonian) to generate the pseudo- 
orbitals. 

To examine the effects of pseudopotentials on the accu- 
racy of the AF QMC calculations, we employ two pseu- 
dopotentials in this study: the optimized LDA-based pseu- 
dopotential^ generated using the OPIUM- 9 package, and the 
HF-based effective core potential developed by Ovcharenko, 
Aspuru-Guzik, and Leste r 30 ' 31 . We will subsequently refer to 
these pseudopotentials as OPIUM and OAL, respectively. The 
semilocal OAL pseudopotentials were converted to the fully 
nonlocal KB form, using the atomic LDA ground state wave 
function as the reference state. The OAL pseudopotential is 
not used in molecules, because it lacks a d projector. An illus- 
tration of this point is given in Table IVlil 

TableUgives parameters describing the OPIUM pseudopo- 
tentials. The second column shows the cutoff energy for each 
atomic species. The same cutoff energies are also used in our 
calculations with OAL pseudopotentials. (The parameters of 
OAL pseudopotentials have been published in Ref.[3(l) These 
were tested for convergence with LDA and then verified with 
AF QMC calculations. 



m. CONVERGENCE STUDIES 

To achieve high accuracy and to minimize the computa- 
tional cost, one should optimize the calculations with respect 
to the number of basis functions, the supercell size, and the 
magnitude of the Trotter time step. In this section, we illus- 
trate the convergence of our method with respect to these pa- 
rameters. 

Ionization energies are defined as IP = E(N — 1) — E(N) 
and IIP = E(N - 2) - E(N), for the singly- and doubly- 
ionized atoms, respectively, where N is the number of elec- 
trons in the neutral atom. The dissociation energy D e is cal- 
culated as the difference between the total energy of the dimer 
at the experimental equilibrium bond length and the energy of 
the isolated atoms, D e = 2E atom — -Edimer- 



A. Planewave convergence 

Convergence with respect to the planewave cutoff energy 
E cul depends on both ffW and in Eq. (fTTat . The ffW 
dependence is similar to that in independent-electron calcu- 
lations. Convergence requires that E cut is sufficient for the 
"hardness" of the pseudopotential and the electronic density 
variations. The dependence has to do with the scattering 
matrix elements in the two-body interaction. In the uniform 
electron gas, for example, requires an E cut given by the 
Fermi energy E-p (for restricted HF), while H ^ will lead to 
a finite convergence error for any finite E cul , which decreases 
as E cut is increased and, for a fixed E cut /Ep, becomes more 
pronounced as the electronic density is decreased. 

Fig. Q] shows the phosphorus atom total energy as a func- 
tion of the planewave cutoff energy E cat for both AF QMC 
and LDA. The calculations were done for a fixed supercell 
size and a pseudopotential whose design cutoff energy is 18 
Ha. In LDA the total energy was converged to within 5 meV 
at this cutoff. The energy decreases monotonically with in- 
creasing E cat in both calculations. We see that the AF QMC 
convergence behavior is similar to LDA, indicating that the 
AF QMC convergence error from is much smaller here 
than that from H W . This trend was found to be typical of the 
systems studied in this paper with the chosen pseudopoten- 
tials. TableQ]shows the cutoff energy for each atomic species. 
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In subsequent calculations for phosphorus, for example, we 
used E cut = 18 Ha in the AF QMC, as indicated by the verti- 
cal dashed line in Fig. Q] 
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FIG. 1: Convergence of the total energy of a phosphorus atom in 
a 14 x 14 x 14 do unit cell. The OPIUM pseudopotential is used 
here. The convergence behavior is similar in both LDA and AF QMC 
methods. A constant shift is added to each data set for convenience, 
such that the converged energies are approximately eV and 0.5 eV 
for LDA and AF QMC, respectively. 



B. Supercell size convergence 

Due to the periodic boundary conditions imposed in the cal- 
culations, the interactions between electrons in the simulation 
cell and their periodic images give rise to finite-size errors. 
To study the behavior of these errors, a series of LDA and 
AF QMC calculations were performed using different system 
sizes for cubic (and some tetragonal) shaped cells. Results 
for the phosphorous atom are shown in Fig. [2] as a function 
of supercell size. The AF QMC energies are seen to converge 
from below while the LDA energy converge from above. This 
is not surprising, since LDA treats the supercell Coulomb in- 
teraction differently from a many-body approach such as AF 
QMC. Fig. [3] shows that the total energy from AF QMC for 
the sulfur atom is nearly a linear function of 1 /Q for this range 
of supercell sizes. At the largest 19 x 19 x 19 dp supercell size, 
the total energy is converged to within about ^0.1 eV. 

To demonstrate the supercell size effect on the energy dif- 
ferences, Fig. |4] shows the calculated dissociation energy of 
P2. The top panel shows the supercell size dependence of 
the dissociation energy, and the bottom panel illustrates the 
convergence error of P and P2 energies for supercells ranging 
from 14 x 14 x 14 eft through 18 x 18 x 18 0%. The total atom 
energy from different supercells deviates no more than 0.2 eV 
from that of the largest supercell. The dimer total energy, on 
the other hand, shows stronger finite-size effect. Most of the 
finite-size error in the dissociation energy thus arises from the 
dimer energy. To avoid any irregular dependence of the en- 
ergy on the aspect ratio (cubic vs. tetragonal supercells), only 
the cubic supercells were used in the extrapolation. 
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FIG. 2: Convergence of the phosphorus total energy with respect to 
the inverse of the simulation cell volume fi. The triangles denote the 
results of LDA calculations while the solid circles denote those of 
the AF QMC. The OPIUM pseudopotential is used here. 
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FIG 3: Convergence of the AF QMC total energy with simulation 
cell size for Sulfur, plotted as a function of the inverse volume, from 
box sizes 12 x 12 x 12 through 19 x 19 x 19 Oq. The OAL 
pseudopotential is used here. The dotted line shows the fitting of the 
total energy as a linear function of 1/fi. The data point shown at 
1/n = shows the extrapolated energy at infinite box size. 



In the ionization energy calculations, the supercells are 
charged +|e| and +|2e| for the X + and X ++ species, re- 
spectively. Charged supercells are ill-defined under peri- 
odic boundary conditions, so an additional neutralizing back- 
ground charge is introduced to maintain charge neutralit y 32,33 . 
As discussed by Makov and Payne 32 ^,, a leading behavior, 
q 2 a/2L, arises from the self-interaction of the neutralizing 
charge with its periodic images, where a is the (supercell- 
dependent) Madelung constant, q is neutralizing charge, and 
L 3 = ft. Correction of the total energy by the leading term 
leads to more rapid size convergence. Figure [5] illustrates this 
effect. The bottom panel shows the slow convergence of the 
total energy with the system size in charged systems, while 
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FIG. 4: Phosphorus dissociation energy (top panel) and the total en- 
ergies of P (times two) and P2 (bottom panel) for different supercell 
sizes, computed using the OPIUM pseudopotential. The dissociation 
energy obtained with the largest simulation cell, 18 x 18 x 18ao is 
within 0.1 eV of the experimental value of 5.08 eV (indicated by the 
horizontal dotted line in the top panel). E com is the energy extrapo- 
lated to 1/fl — > (using cubic cells only). 



the upper panel shows the more rapid convergence after the 
correction has been made, i.e., the slowly convergent q 2 a/2L 
contribution has been subtracted. 



Trotter-Suzuki decomposition in Eq. ( fTol l. The Trotter er- 
ror can be eliminated by extrapolation, as demonstrated in 
Fig. [6] for P, P + , and P ++ . The results reported in the next 
section use either a linear extrapolation a fixed At (~ 0.05 
Ha -1 ) which is sufficiently small so that the Trotter error is 
well within the statistical error. 



-177 
-177 
-177 
-177 

'o -166 
>>-167 

DC 

| -167 



o 
t— 



1 I 


I 


? 






_ 








*-* , 












- p + 








i 


_ 

: p ++ 


— a- 


i 


—5- 


i 

1 J — 



0.05 0.10 

Time step Ax (units of Ha ) 

FIG. 6: Trotter errors for P, P + and P ++ energies (shown in the 
top, middle, and bottom panels, respectively) in a 14 x 14 x 14 clq 
simulation box. The data points at At = represent the Trotter 
extrapolation. 
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FIG. 5: An illustration of the size effect in the calculation of charged 
atoms. The top panel shows the corrected total energies for P, P + and 
P ++ , while the bottom panel shows the uncorrected ones. Cell sizes 
used in this study are the same as in Fig. [4] The errorbars are smaller 
than the point size. 



C. Trotter time step error 

The Trotter error arises from neglecting higher order terms 
of the imaginary-time propagator, e~ , when we apply the 



IV. RESULTS 

Table [TT] shows dissociation energies from planewave AF 
QMC calculations compared with experimental values and 
with LDA, GGA, HF, and diffusion Monte Carlo (DMC) 
calculated results. LDA trial wave functions were used 
in the molecular calculations, corresponding to the fol- 
lowing electronic configurations: HL^SsHL^fp f° r P2, 
HL^HWp^p for S 2< ^SsHWp^p for Cl 2 , and 
^Is^Is ^Ip^lp f° r A s 2- All dimers except S2 have "closed- 
shell" configurations in which all occupied orbitals are fully 
filled. All the calculations used OPIUM pseudopotentials. 

Both restricted and unrestricted trial wave functions were 
tested in the AF QMC calculations, but there were no differ- 
ence within statistical errors. (In restricted trial wave func- 
tions, the orbitals with minority spin are identical to the corre- 
sponding majority-spin orbitals.) Unrestricted HF trial wave 
functions were also used to calculate the dissociation energy 
of S2, which has an open-shell configuration. As shown in 
Table HIl there is no difference within statistical errors. 

The overall agreement between the AF QMC results and 
experiment is very good. The LDA and GGA slightly over- 
estimate the dissociation energy, while the HF method signif- 
icantly underestimates it. The heaviest dimer we calculated, 
AS2, is also in excellent agreement with experiment and com- 
pares favorably with results from other quantum chemistry 
methods^. 

We show in Table [III] the first and second ionization en- 
ergies of Al, Si, P, S and CI. For comparison, experimen- 



8 



TABLE II: Calculated dissociation energies (in eV) using LDA, AF 
QMC, DMC, HF, and GGA methods. Experimental values are in the 
last column (with the zero-point energy removed). The AF QMC cal- 
culations used LDA trial wave functions, except for sulfur, where we 
have also used unrestricted HF trial wave functions. Values from 
the largest supercell are shown here, since convergence has been 
reached. The statistical errors are given in parentheses. The HF re- 
sults are obtained from an in-house planewave based code. The AF 
QMC result for Si2 is taken from Ref.0. 

LDA AF QMC (LDA) DMC 34 HF GGA 35 Expt^ 



Si 2 
P> 

S, 

Cb 
As 2 



3.88 
5.97 
5.61 

3.12 
5.04 



3.12(8) 

5.19(16) 

4.63(17) 

4.48(19)° 

2.78(10) 

3.97(17) 



4.73(1) 
4.31(1) 



1.74 
2.29 



3.21(13) 
5.22 5.08 
4.94 4.41 



2.38(1) 0.67 2.76 



2.51 
3.96 



"Using HF trial wave function. 



tal values and results from LDA and HF calculations are also 
shown. The AF QMC calculations are performed using both 
the OPIUM and OAL pseudopotentials together with LDA 
trial wave functions, except for sulfur with the OAL pseu- 
dopotential, where HF trial wave functions are also used. 
Again, no statistically significant dependence on the trial wave 
function is seen. The tabulated results are converged with 
respect to size effects, which are negligible compared to the 
QMC statistical error for supercells larger than 16x16x16 Oq. 

Ionization energies obtained using the LDA are generally in 
very good agreement with experiment, and the present results 
are consistent with this, with deviations typically within 0.1 
eV for IP and 0.3 eV for IIP, regardless of which pseudopo- 
tential is used. The HF method results in larger deviations of 
0.5-1.3 eV compared to the experiment. 

AF QMC results for Al and Si are in good agreement with 
experiment, and for Si the agreement is comparable to that 
obtained using DMC 4 !, 8.166(14) eV and 24.444(22) eV for 
IP and IIP, respectively. For P, S, and CI, however, the agree- 
ment between AF QMC and experiment is not as uniform. In 
particular, there is a significant dependence on the choice of 
the pseudopotential. For P and S, the AF QMC ionization en- 
ergies are better estimated using OAL, while for CI, OPIUM 
pseudopotential gives better results. We will discuss this de- 
pendence in the next section. Agreement between the best AF 
QMC and experiment values is in general very good. 



V. DISCUSSION 

For dissociation energies, the agreement between AF QMC 
and experiment was uniformly very good. The appearance 
of larger discrepancies between AF QMC and experiment for 
ionization energies is somewhat surprising. The AF QMC cal- 
culations above were systematically converged with respect to 
finite size effects, Trotter time step error, and planewave basis 
size. Two remaining possibilities, errors arising from the use 
of the phaseless approximation and errors due to the use of 
pseudopotentials, are discussed in this section. 



Our overall experience with the phaseless AF 
QM C 4 ' 9,10,12 ' 24 suggests that the error due to the phase- 
less approximation itself is typically small. Recent AF 
QMC calculations using a Gaussian basis^S show that, in 
a variety of atoms and molecules, the AF QMC agrees 
well with experiment or high-level quantum chemistry 
methods such as the coupled cluster with single and double 
excitations and perturbative corrections for triple excitations 
[CCSD(T)] 42 i 43 i 44 - 45 . Open-shell systems such as P+, P++, 
S, S ++ , CI, and Cl + , where the p shell is neither half- nor 
fully-filled, tend to be more difficult to treat in general 4 ^. In 
these cases a single-determinant trial wave function breaks 
the symmetry by using only one of the degenerate states, 
and the phaseless approximation could affect the accuracy 
of the results. Some indication of this may be present in 
our results (see Table HIB where larger errors are observed 
for energy differences between half-filled and "open-shell" 
systems, such as the P + ionization energy. It is possible to use 
multideterminant trial wave functions in these cases. We have 
done several tests, in which we "symmetrize" the trial wave 
function, resulting in a linear combinations of three different 
determinants. This is designed to equally treat the p x , p y , 
and p z orbitals in the open shell. However, there seems to be 
no observable improvement over the single-determinant trial 
wave function at the level of statistical accuracy in this paper. 

To isolate and quantitatively evaluate the errors due to 
the phaseless approximation on the second-row atoms stud- 
ied here, we performed calculations with both AF QMC and 
CCSD(T) methods for CI, C1+, and C1++, using identical 
Gaussian basis sets and the OAL pseudopotential. (The OAL 
pseudopotential was chosen since its form is already com- 
patible with standard quantum chemistry programs.) Thus, 
both the AF QMC and CCSD(T) methods are applied to the 
same many-body Hamiltonian, expressed in the Hilbert space 
spanned by the selected Gaussian basis set. The CCSD(T) 
method is approximate but is known to be very accurate in 
atoms and in molecules near equilibrium. For the compar- 
ison, we employed an uncontracted aug-cc-pVDZ 4 -^ basis 
set, where Gaussian functions with exponents larger than 98 
are removed, resulting in a (7s7p2d) basis set. As shown in 
Table HVl the AF QMC and CCSD(T) absolute total energies 
agree to within 0.07 eV, and their ionization energies agree to 
within 0.04 eV. These results indicate that the intrinsic error 
due to the phaseless approximation in the above planewave 
calculations is likely also small, consistent with previous ex- 
perience with phaseless AF QM C 4 ' 9 ' 10 ' 12 ' 24 . This suggests that 
errors due to the use of pseudopotentials are largely responsi- 
ble for the deviations in ionization energies noted above. 

The pseudopotentials used here were generated using 
independent-electron HF or DFT mean-field type calculations 
of atomic reference systems. While the transferability of these 
pseudopotentials to HF or DFT calculations of molecules or 
solids is well understood, their accuracy in many-body calcu- 
lations is more problematic i 4 -^ The dependence of the AF 
QMC results in Table IllTl on the choice of pseudopotentials is 
consistent with this. 

To estimate the pseudopotential errors in our AF QMC cal- 
culations and to obtain insight into their origin, we have car- 
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TABLE III: First and second ionization energies for several atoms, in eV. We employ two pseudopotentials, the LDA-based OPIUM pseudopo- 
tential and HF-based effective core potential (OAL). In most cases, we use LDA trial wave functions (see text). The LDA calculations were 
done using ABINIT. The HF results were obtained using an in-house HF program for OPIUM pseudopotential and GAUSSIA N 39 ' 40 for OAL. 
Experimental values are taken from Ref.l36l The quantity AIP below is the difference between the AF QMC-calculated and the experimental 
ionization energies, for which we show the average and r.m.s. average over all the species shown here. 







IP(X 








IIP (X 


-» A++) 




Pseudo 


HF 


LDA 


AF QMC 


Expt 


HF 


LDA 


AF QMC 


Expt 


Al OAL 




5.88 


5.88(2) 


5.99 




24.46 


24.66(2) 


24.81 


Si OPIUM 






8.18(2) 


8.15 






24.59(4) 


24.50 


P OPIUM 


9.97 


10.57 


10.74(6) 


10.49 


29.41 


30.42 


30.79(6) 


30.26 


OAL 


9.94 


10.41 


10.61(3) 




29.11 


29.97 


30.28(6) 




S OPIUM 


9.33 


10.45 


10.09(7) 


10.36 


32.42 


33.85 


34.16(7) 


33.67 


OAL 


9.21 


10.35 


10.08(2) 




32.06 


33.51 


33.61(2) 




CI OPIUM 


11.69 


13.12 


12.96(11) 


12.97 


34.36 


36.92 


36.76(10) 


36.78 


OAL 


11.76 


13.02 


12.89(6) 




34.30 


36.64 


36.25(6) 




(AIP(OPIUM)) 






0.00 








0.27 




(AIP(OPIUM)) 


rms 




0.19 








0.36 




(AIP(OAL)) 






-0.09 








-0.18 




<AIP(OAL)) rms 






0.16 








0.28 





TABLE IV: Calculated total energies (E) and ionization energies (IP) 
for CI, using AF QMC with a Gaussian basis together with the cor- 
responding CCSD(T) results, using identical basis sets. The OAL 
pseudopotential is used with a double-zeta quality (7s7p2d) Gaus- 
sian basis set. All energies are in eV. Since CCSD(T) is known to 
be accurate for atoms, differences compared to Gaussian-based AF 
QMC provide an estimate of errors due to the phaseless approxima- 
tion in AF QMC. 





AF QMC 


CCSD(T) 




E IP 


E IP 


CI 


-403.91(1) 


-403.96 


C1+ 


-391.37(1) 12.54(2) 


-391.44 12.51 


C1+ 4 


-368.43(1) 35.49(2) 


-368.43 35.53 



ried out several additional calculations. We first performed 
pseudopotential and all-electron (AE) CCSD(T) calculations 
to estimate the transferability of the pseudopotential in a 
many-body context. First and second ionization energies for 
P, S, and CI atoms were calculated using both methods. The 
coupled-cluster calculations were performed using the GAUS- 
SIAN98 and GAUSSIAN03 packages^il£. To eliminate basis set 
convergence errors, we performed a series of AE calculations 
using the aug-cc-pwCVxZ basis set o 47 ' 48 ' 50 , where x — D, T, 
Q for double, triple, and quadruple zeta basis sets, respec- 
tively. The infinite-basis estimate of the total energy, E^, is 
then obtained using extrapolation^ 

E^^Ex- be- cx , (32) 

where x is 2, 3, 4 for double, triple, and quadruple zeta basis 
sets, respectively, and b and c are fitting parameters. We then 
take the difference of the extrapolated energies as the ioniza- 
tion potential shown in Table [Vl For the pseudopotential cal- 
culations, the OAL ECP is used. Here we use the aug-cc- 
pWxZ basis set a 47 ' 48 that are fully uncontracted, and again we 



use the extrapolation scheme in Eq. d32l . 

TABLE V: Estimates of pseudopotential errors: CCSD(T) results for 
the first and second ionization energies of P, S, and CI. All energies 
are in eV. OAL and AE results are shown together with the error in 
the ionization energy due to the pseudopotential, AIP psp . 







IP(X 








IIP(X 








OAL 


AE 


AIP pS p 


Expt. 


OAL 


AE 


AIIP psp 


Expt. 


p 


10.48 


10.47 


0.01 


10.49 


30.13 


30.18 


-0.05 


30.26 


S 


10.24 


10.29 


-0.05 


10.36 


33.64 


33.64 


0.00 


33.67 


CI 


12.92 


13.05 


-0.13 


12.97 


36.51 


36.77 


-0.26 


36.78 



The CCSD(T) results for AE and OAL calculations are pre- 
sented in Table [V] AE CCSD(T) ionization energies are seen 
to be in excellent agreement with experimental values. [We 
note that AE calculations using aug-cc-pCVxZ basis set (a 
variation of aug-cc-pwCVxZ) yield almost identical result, 
where the estimated IP differs by < 0.03 eV, and IIP by 
< 0.07 eV.] The difference between the AE and OAL results 
provides an estimate of the error due to the pseudopotential in 
many -body calculations. The results in Table Pvl indicate that 
the OAL pseudopotential tends to underestimate the ioniza- 
tion energies. The rms error attributable to the pseudopoten- 
tial is about 0.12 eV. This is not negligible compared to the 
overall AF QMC error, which is (AIP(OAL)) rms = 0.25 eV 
for P, S, and CI first and second ionization energies. 

An additional possible source of pseudopotential error in 
the planewave AF QMC calculations is the use of the fully 
nonlocal, separable Kleinman-Bylander (KB) construction of 
the pseudopotential. For example, the OAL ECP is defined in 
the usual semilocal form, which is used in quantum chemistry 
programs: 

Vt\r) = Y,\ Y irn)V l (r)(Y lm \. (33) 

l,m 
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where Vi(r) is the angular-momentum-dependent potential. 
For efficient use in the planewave calculations, it is common 
to express this pseudopotential in the fully nonlocal separa- 
ble KB form shown in Eq. (13 lab . While the KB and semilo- 
cal forms are identical when they act on the reference atomic 
state, the KB form can differ for other states. 

To investigate the effect of pseudopotential KB formation, 
we perform LDA^ 2 calculations with the OAL pseudopoten- 
tial: (1) planewave basis calculations with the KB form of the 
OAL using the ABINIT package (calculations are converged 
with respect to E cut and box size), and (2) local basis calcula- 
tions with the semilocal form of the OAL using GAUSSIAN98 
(again we use the sequence of aug-cc-pV.xZ basis sets to ex- 
trapolate to the infinite basis limit). The OAL pseudopotential 
was converted to the KB form using pseudo-orbitals obtained 
in an LDA calculation for the neutral atom. For the purpose 
of constructing the fully nonlocal projectors, the effects of us- 
ing LDA rather than HF pseudo-orbitals are not expected to 
to be significant. In both methods, the total-energies are con- 
verged to within 0.5 mHa (« 0.01 eV). Table lVII presents the 
results for the CI ionization energies. The OAL pseudopoten- 
tial expressed in the KB form tends to underestimate the LDA 
total energies, with the discrepancy increasing with the ion- 
ization state. Ionization energies are thus underestimated by 
up to 0.15 eV. The same trend is observed in the calculated 
planewave AF QMC ionization energies compared to experi- 
ment, which indicates that the KB form may contribute errors 
of the order of 0.1 — 0.2 eV for CI using the OAL pseudopo- 
tential. It is clear that the quality of the pseudopotential is 
crucial for obtaining accurate results with AF QMC. 

TABLE VI: LDA calculations of chlorine energies, the Kleinman- 
Bylander (KB) and semilocal (SL) forms of the OAL pseudopoten- 
tial. The programs used are indicated in parantheses, and all energies 
are in eV. 





Kleinman-Bylander 


Semilocal 






(ABINIT) 


(GAUSSIAN98) 


(E KB -E SL ) 


System 


E IP 


E IP 


CI 


-403.798 


-403.752 


-0.05 


C1+ 


-390.842 12.956 


-390.701 13.051 


-0.14 


C1++ 


-367.217 36.581 


-367.021 36.731 


-0.20 



These test calculations are limited in scope. Small sys- 
tematic errors may well be still present, for example from 
the planewave size-extrapolations, from Gaussian basis set 
extrapolations, and from approximations inherent in coupled 
cluster calculations at the CCSD(T) level. They suggest, how- 
ever, a rather consistent picture for understanding the AF 
QMC results in TablellJTl Pseudopotential errors due to differ- 
ent origins appear to be the main cause for the discrepancies 
with experimental values. When such errors are removed, it 
seems that the accuracy of the planewave AF QMC is at the 
level of 0.1 eV. 

Since the quality of a pseudopotential within HF and LDA 
can be easily determined (by its ability to reproduce AE re- 
sults), it is interesting to see how well this correlates with the 
performance of the pseudopotential in a many-body calcula- 
tion. In Table I VIII we show various energies obtained from 



TABLE VII: Phosphorus ionization and dissociation energies for P , 
P ++ , and P 2 computed using LDA, HF, AF QMC, and CCSD(T) 
methods, shown below in eV. AE results are provided to benchmark 
the pseudopotentials. The numbers shown in boldface are those clos- 
est to the AE results (for LDA and HF) or the experimental values 
(for AF QMC). 



PSP 


LDA 


HF 


AF QMC 


CCSD(T) 


IP(P^ 


P + ) 






expt = 10.49 


OPIUM 


10.57 


9.97 


10.74(6) 




OAL 


10.41 


9.94 


10.61(3) 


10.48 


AE 


10.53 


9.91 




10.47 


IIP (P - 


-P++) 






expt = 30.26 


OPIUM 


30.42 


29.41 


30.79(6) 




OAL 


29.97 


29.11 


30.28(6) 


30.13 


AE 


30.37 


29.08 




30.18 


D e (P 2 - 


-> 2P) 






expt = 5.08 


OPIUM 


5.97 


1.74 


5.19(16) 




OAL 


5.29 


0.98 


3.88(8) 


4.39 


AE 


6.18 


1.65 




4.98" 



"Frozen-core calculation. 



LDA, HF, AF QMC, and CCSD(T) calculations in phospho- 
rus. (Similar trends also hold for sulfur and chlorine.) Re- 
sults are shown using the OPIUM and OAL pseudopoten- 
tials, together with AE results for independent-electron and 
CCSD(T) methods. In the LDA calculations, the LDA-based 
OPIUM potential performs uniformly better for all quantities. 
In the HF calculations, the HF-based OAL pseudopotential 
performs well, except for the dissociation energy. (As men- 
tioned earlier, the OAL pseudopotential lacks a d projector 
in its construction, which results in poor molecular energies.) 
Overall, the HF method appears to give a better indication of 
pseudopotential performance with AF QMC. In the chlorine 
second ionization energy, for example, HF predicts that the 
OPIUM pseudopotential performs better than the OAL pseu- 
dopotential: IIPhf, opium = 34.37 eV, IIPhf, oal = 34.30 eV, 
compared to IIPhf, ail-electron = 34.36 eV. AF QMC results in 
Table Hn] show a similar trend. We also note in Table [TTT1 that 
the LDA-based OPIUM pseudopotentials tend to overestimate 
the ionization energies, while the OAL pseudopotentials do 
the opposite. The tabulated rms averages suggest that the 
performance of AF QMC with the LDA-generated OPIUM 
pseudopotential varies more widely across different atomic 
species, especially in the second ionization energies. The HF- 
generated OAL ECP, on the other hand, performs more consis- 
tently and yields better agreement with experiment in the ma- 
jority of species studied here. Testing pseudopotentials using 
Hartree-Fock calculations may, therefore, be a useful predic- 
tor of their performance in the many-body AF QMC method. 



VI. SUMMARY 

We have presented electronic structure calculations in 
atoms and molecules using the phaseless AF QMC method 
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with a planewave basis and norm-conserving pseudopoten- 
tials. Various algorithmic issues and characteristics were de- 
scribed and discussed in some detail, and we have illustrated 
how the AF QMC method can be implemented by utilizing 
standard DFT planewave techniques. The structure of the 
AF QMC calculation is an independent collection of random 
walker streams. Each stream resembles an LDA calculation, 
which makes the overall computational scaling of the method 
similar to LDA calculation with a large prefactor. This makes 
the AF QMC approach more efficient than explicit many- 
body methods. All of the reported results were obtained us- 
ing single-determinant trial wave functions, directly obtained 
from either LDA or HF calculations. This reduces the demand 
for wave function optimization in QMC and is potentially an 
advantage. The method also offers a different route to the sign 
problem by carrying out the random walks in Slater determi- 
nant space. Because our method is based in Slater determinant 
space, any single-particle basis can be used. 

Results for the dissociation and ionization energies of 
second-row atoms and dimers in Group 3A-7A, Al, Si, P, 
S, CI, as well as the As 2 dimer, were presented using the 
planewave-based phaseless AF QMC method. The effects of 
the phaseless approximation in AF QMC were studied, and 
the accuracy of the pseudopotentials were examined. Com- 



parisons were made with experiment and with results from 
other methods including the DMC and CCSD(T). Errors due 
to the phaseless approximation were found to be small, but 
non-negligible pseudopotential errors were observed in some 
cases. In addition to pseudopotential errors that arise due to 
their construction in mean-field type DFT or HF calculations, 
possible errors in ionization energies arising from the sepa- 
rable Kleinman-Bylander form of the pseudopotentials were 
also observed. With the appropriate pseudopotentials, the 
method yielded consistently accurate results. 
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